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Preface 


This  report  is  the  result  of  my  investigation  of  the 
Explicit,  Pure  Implicit,  Crank-Nicolson  and  Douglas  finite- 
Difference  methods  for  solution  of  the  one-dimensional 
transient  heat -conduction  equation  in  inhomogeneous  material 

Initially,  I  sought  to  solve  three  variations  of  the 
basic  problem  for  several  values  of  the  Fourier  Modulus, 

At/ (Ax)2,  and  test  two  or  three  methods  for  handling  the 
inherent  discontinuity  found  in  the  problem.  However, 
time  and  computer  assets  allowed  for  complete  solution  of 
only  two  variations  of  the  problem  for  four  values  of  the 
Fourier  modulus,  with  one  treatment  of  the  discontinuity. 
Yet,  the  value  of  traditional  finite-difference  schemes 
over  newer,  more  sophisticated  schemes  for  engineering  work 
was  still  substantiated. 

I  would  like  to  express  my  appreciation  to  Dr.  Bernard 
Kaplan  of  the  Air  Force  Institute  of  Technology  for  his 
guidance,  and  to  Dr.  W.  Kessler  of  the  Air  Force  Materials 
Laboratory  for  sponsoring  this  research  project.  (Also,  to 
Sharon  Gabriel  for  doing  such  a  magnificent  and  professional 
job  with  the  typing  that  I  can't  miss  getting  a  passing 
grade! ) 


Kenneth  W.  Blevins 
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Abstract 

The  transient  heat  conduction  equation  is  solved  for 
inhomogeneous  media  using  the  Explicit,  Pure-Implicit, 

Crank -Nicolson  and  Douglas  finite-difference  methods,  and 
the  numerical  solutions  are  investigated  with  respect  to 
accuracy  and  stability.  The  inherent  discontinuity  between 
the  initial  and  boundary  conditions  is  accounted  for  by  mesh 
refinement.  For  the  two  versions  of  the  problem  for  which 
the  four  numerical  methods  are  investigated,  all  four  methods 
are  found  to  be  of  equivalent  accuracy  for  small  values  of 
the  Fourier  Modulus,  At/ (Ax)2~^>>  While  the  Pure- Implicit , 
Crank-Nicolson  and  Douglas  methods  are  unconditionally  stable, 
the  Crank-Nicolson  and  Douglas  methods  are  very  inaccurate 
at  large  values  of  the  Fourier  Modulus  due  to  oscillatory 
behavior . 


A  COMPARISON  OF  FINITE-DIFFERENCE  METHODS 
FOR  THE  SOLUTION  OF  THE  TRANSIENT  HEAT 
CONDUCTION  EQUATION  IN  INHOMOGENEOUS  MEDIA 

I.  Introduction 


General 

Background .  Very  few  practical  engineering  problems 
involving  partial  differential  equations  can  be  solved  in 
closed  form.  The  difficulty  may  be  due  to  irregular  or  mixed 
geometry,  or  complicated  boundary  conditions.  Further, 
closed  form  solutions  that  are  obtainable  usually  contain 
infinite  series,  special  functions,  or  transcendental 
equations  for  eigenvalues  such  that  the  numerical  evaluation 
of  an  analytic  solution  itself  may  be  a  difficult  task 
(Ref  16:3).  Hence,  numerical  approximation  methods  are 
used  which  provide  adequate  numerical  solutions  simply  and 
efficiently.  Finite  differences  have  been  in  use  at  least 
since  1768  (Ref  1:1),  and  are  the  most  frequently  used  and 
universally  applicable  numerical  approximation  method  for 
solving  partial  differential  equations  (Ref  20:4). 

Various  finite  difference  schemes  of  different  degrees 
of  sophistication  are  given  in  the  literature  for  the 
numerical  solution  of  the  transient  heat -transfer  equation. 

For  simplicity  of  presentation  and  to  maintain  generality, 
these  schemes  are  usually  presented  for  homogeneous  materials. 
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Much  work  has  been  performed  comparing  the  accuracy  and 
efficiency  of  these  different  schemes  in  solving  the  homo¬ 
geneous  material  transient  heat -transfer  problem.  However, 
the  solution  of  many  important  engineering  heat  transfer 
problems,  such  as  those  involving  re-entry  bodies  and  jet 
engine  nozzles,  requires  handling  non-homogeneous  (composite) 
materials.  Some  of  the  finite-difference  schemes  available 
for  use  may  be  potentially  more  advantageous  with  respect 
to  accuracy  and  computer  time  than  those  currently  employed 
in  engineering  practice. 

Finite  Difference  Formulation 

When  using  the  finite  difference  technique  to  solve  a 
partial  differential  equation  (PDE) ,  a  network  of  grid  points 
must  first  be  established  throughout  the  region  of  interest. 
The  grid  lines  are  drawn  parallel  to  the  space  and  time  axis 
and,  for  one-dimensional  spatial  coordinates,  are  separated 
by  Ax  and  At  .  The  intersections  of  the  grid  lines  are 
the  grid  points  or  nodes  and  are  identified  by  spatial  and 
time  coordinates.  Let  x  and  t  be  the  independent  varia¬ 
bles  of  the  PDE  to  be  solved,  and  let  the  grid  spacings  be 

Ax  for  space  and  At  for  time.  The  spatial  coordinate  of 

an  arbitrary  node  will  be  B  and  the  spatial  coordinates  of 

the  adjacent  nodes,  A  and  C  .  The  value  of  U  at  a 
particular  node  will  be  specified,  at  time  t  ,  by  a  sub¬ 
script  and  at  time  t+6t  by  the  same  subscript  and  a  prime 
superscript,  as  depicted  in  Figure  1. 
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Figure  1.  Grid  Points  and  Notation 

Assume  the  exact  solution  to  the  PDE  is  UA  s  =  U(x,t) 
and  let  the  approximation  to  the  exact  solution  be  up  p  g  ‘ 
Next,  suitable  finite  difference  expressions  must  be  found 
which  will  approximate  the  partial  derivatives  of  the  original 
PDE.  After  substituting  the  finite  difference  expressions 
into  the  original  PDE,  the  resulting  finite-difference  equation 
(FDE)  is  solved.  However,  the  solution  of  the  FDE  is  only  an 
approximation  to  the  actual  solution  of  the  original  PDE 
because  the  derivatives  have  just  been  approximated  by  dif¬ 
ference  quotients  over  small  intervals.  As  the  size  of  the 
intervals  approaches  zero,  the  difference  quotients  approach 
the  value  of  the  derivatives  they  approximate. 

The  procedure  used  in  this  study  for  deriving  the  finite 
difference  equations  consists  of  approximating  the  derivatives 
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of  the  PDE  by  a  truncated  Taylor  series.  To  do  this,  first 
refer  to  the  grid  points  shown  in  Figure  1.  For  a  node  U„  , 

D 

located  halfway  between  nodes  and  U^,  such  that 

Ax  =  Xg  -  -  Xg  and  At  =  t'  -  t  ,  the  Taylor  series 

expansion  about  Ug  along  the  x  coordinate  system  gives 


UA  = 

dUB  1 

UB  "  Ax  "35T  +  7 

(AX)2 

d!UB 

dx2 

(1) 

UC  " 

dUB  1 

(AX)2 

d’UB 

dx2 

(2) 

where  terms  of  order  3,  0(Ax3)  ,  have  been  neglected.  Adding 

and  subtracting  the  two  equations  and  rearranging  terms 


dUB  .  "C  '  UA 
dx  2  (Ax) 


(3) 


d2UB  .  UC  -  2UB  *  UA 
dx2  (Ax)2 


(4) 


Equations  (3)  and  (4)  are  respectively  referred  to  as  first 
and  second  central  difference  equations.  Similarly,  where 
lT(x,t)  =  U(x,t  +  At)  ,  the  first  forward  difference  equation 
for  the  time  variable  is 


dUB  UB  '  UB 
ITT  =  — At — 


(5) 
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Approximating  a  PDE.  The  partial  differential  equation 


3U  =  3^U 

8t  =  3x2 


(6) 


may  be  approximated  by  substituting  (4)  and  (5)  into  (6) 
to  get 


UB  -  UB 
At 


UA  -  2UB  +  UC 
(Ax)2 


(7) 


which  is  a  formula  for  the  unknown  value  U'  in  terms  of 

15 

the  known  values  ,  Ug  ,  Uc  .  Referring  to  Figure  2a, 

solving  for  one  unknown  directly  in  terms  of  known  values 
is  called  an  explicit  method  (Ref  20:13),  and  (7)  is  the 
Euler  or  Explicit  formula. 

A  different  approximation  to  (6)  is  given  by 


UB  -  UB 
At 


-  2ub  *  u6 

(Ax)2 


(8) 


which  is  generally  known  as  the  Pure-Implicit  formula. 
Referring  to  Figure  2b,  ,  Ug  ,  are  unknown  and 

must  be  solved  for  in  terms  of  Ug  which  is  known.  This 
leads  to  a  set  of  n  -  1  equations  in  n  -  1  unknowns 
where  n  +  1  is  the  number  of  grid  points  and  the  solution 
is  known  at  the  boundaries.  The  equations  may  be  solved  by 
straightforward  Gaussian-elimination.  However,  when  the 
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Figure  2.  Schematic  Representation  of  Explicit 
and  Implicit  Numerical  Methods 

matrix  of  the  coefficients  is  written,  the  non-zero  terms 
align  themselves  along  the  three  main  diagonals  of  the  matrix, 
forming  a  tri-diagonal  matrix,  and  a  special  adaptation  of 
the  Gaussian-elimination  procedure,  the  Thomas  method,  may 
be  used  for  the  solution  (Ref  18:199). 

All  forms  of  finite  difference  equations  may  be  classi¬ 
fied  as  explicit  or  implicit.  Explicit  formulae  march  forward 
in  time  as  the  solution  at  each  present  point  is  found  in 
terms  of  preceding  values,  while  implicit  formulae  involve 
the  solution  of  simultaneous  equations  of  present  unknown 
values  in  terms  of  preceding  known  values  (Ref  1:42). 

General  Finite  Difference  Equation.  A  general  finite- 
difference  form  of  Eq  (6)  has  been  given  by  Crandall  (Ref  6: 
318)  as  a  weighted  combination  of  the  explicit  and  implicit 
relations  pictured  in  Figure  2,  and  is  given  by 
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M(U'  -  UB) 


-  £(u;  -  2Ug  <■  U')  *  (l-f)tUA  -  2Ub  *  Uc) 


where 

U  * 

the  approximate  temperature  at 

point  as  shown  in  Figure  1 

a  specified  grid 

M  « 

the  inverse  Fourier  Modulus,  a 

ratio  (Ax)2/ (At) 

dimensionless 

f  = 

weighting  factor. 

In  general  practice,  0  s  f  s  1.0  .  A  schematic  for  the 
general  equation  is  given  in  Figure  3. 


Figure  3.  Schematic  Representation  of  the 

General  Finite-Difference  Equation 

Table  1  includes  several  values  of  f  and  common  nomencla¬ 
ture  for  the  finite-difference  schemes  thus  determined. 


(9) 
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Table  1 

Finite  Difference  Schemes  Determined 
by  Values  of  the  Weight  Factor  f 


Value  of  f 

Finite  Difference  Scheme 

0 

Explicit 

1/2 

Crank-Nicolson 

1 

Pure- Implicit 

lf,  (Ax)  2  . 

6a(x)At' 

Douglas 

Obj  ective 

The  basic  purpose  of  this  thesis  is  to  investigate  the 
relative  accuracy  of  the  Explicit,  Pure  Implicit,  Crank- 
Nicolson,  and  Douglas  finite  difference  schemes  used  to  rep¬ 
resent  the  parabolic  partial  differential  equation  for 
transient  heat  conduction  in  non-homogeneous  media,  and 
to  compare  the  methods  for  ease  and  efficiency  of  computation. 

Approach  and  Scope 

A  straightforward  approach  to  investigating  the  accuracy 
of  numerical  methods  is  to  solve  a  given  problem  analytically 
and  numerically,  and  compare  the  errors  of  the  methods  used. 

This  is  the  general  approach  used  in  this  thesis.  Although 
it  would  be  ideal  to  apply  these  methods  to  a  general  two- 

I 

< 

i 
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dimensional  or  three-dimensional  time-dependent  problem,  and 
compare  the  solutions  obtained  to  the  analytic  solution  of 
the  same  problem,  the  analytic  solution  to  such  a  problem  is 
extremely  difficult,  if  at  all  possible  to  obtain.  Thus,  this 
investigation  was  restricted  to  the  study  of  a  one-dimensional 
problem  as  found  in  Carslaw  and  Jaeger  (Ref  4:312-313).  Due 
to  the  limited  usefulness  and  severe  stability  restrictions 
of  the  explicit  form,  as  shown  in  Chapter  II,  the  main  thrust 
of  my  effort  was  to  compare  the  Crank-Nicolson,  Douglas  and 
Pure  Implicit  forms  of  solution. 

Computer  Program 

Investigation  of  the  discussed  numerical  methods  and 
numerical  computation  of  the  analytic  solution  was  done  with 
the  use  of  computer  programs  written  by  the  author  in  FORTRAN 
Version  5.  The  programs  were  executed  on  the  Aeronautical 
Systems  Division  computer  at  Wright-Patterson  Air  Force  Base, 
Ohio. 
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II.  Theor 


The  Physical  Problem 

The  heat  transfer  problem  investigated  is  illustrated 
in  Figure  4. 


The  problem  consists  of  a  solid  which  extends  to  infinity 
in  the  x-direction.  The  material  contains  no  heat  sources. 
The  entire  solid  is  initially  cooled  to  0°.  Subsequently, 
the  temperature  at  the  surface  of  the  solid  at  x  equal  to 
zero  is  immediately  raised  to  and  maintained  at  a  temperature 
of  one  degree.  The  temperatures  used  have  no  significance 
other  than  to  simplify  the  initial  and  boundary  conditions. 
All  other  faces  of  the  solid  are  insulated,  thus  heat  is 
transferred  only  in  the  x-direction.  The  conductivity  of 
the  material  is  proportional  to  x  to  the  nth  power,  where 


n  is  greater  than  zero  and  less  than  unity,  (0  <  n  <  1)  . 
The  material  used  has  no  bearing  on  the  investigation,  and 
the  physical  parameters  of  specific  heat  and  density  have 
been  considered  constant.  The  governing  partial  differential 
equation  for  this  conduction  problem  is  given  by 

(IC  5?)  -  pc  at  0 


where 

U 

t 

x 

P 

c 

K 


the  exact  temperature  at  a  specific 

node  and  time 

time 

node  location 
density 
specific  heat 

conductivity  at  a  specific  node  such  that 


K  =  KQxn  0  <  n  <  1 

K  =  constant 
o 

The  initial  and  boundary  conditions  which  must  be  satisfied 
are 

U  (x ,0)  =  0 

U  (0,t)  =  i 

U  (»,t)  =  0  (11) 
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Exact  Solution 


The  exact  solution  to  this  problem  is  given  by  Carslaw 
and  Jaeger  (Ref  4:412-413)  as 

°o 

U(x»t)  =  nvr  /  e"y  ^v"1  ^  uz) 

X 


where 

U  =  temperature  at  node  x  and  time  t 
y  =  variable  of  integration 
v  =  (1  -  n)/(2  -  n)  0  <  n  <  1 

T(v)  =  Gamma  function  evaluated  at  v 
X  =  (x2'n)/((2  -  n) 2kt) 

k  =  KQ/pc  where  Kq,  p,  c  were  defined  for  (10). 

The  integral  in  Eq  (12)  must  be  evaluated  numerically.  The 
approximation  used  for  this  computation  was  given  by 
Gradshteyn  and  Rhyzik  (Ref  9:940-941)  as 


e-1  t"'1  dt 


r(oO  -  l 

n=0 


(-l)n  x^*1^ 


Figure  5  depicts  the  exact  analytical  solution  for  four 
variations  of  the  problem  where  n  =  0.0  ,  0.25  ,  0.5  ,  and 
0.75  ,  and  Figure  6  depicts  the  temperature  profile  at  several 
times  for  n  =  .25  . 
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Figure  5.  Exact  Solution  for  n=0, 

.25,  .5,  and  .75  at  x=.05 


TEMPERATURE  PROFILES 
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X 

Figure  6.  Temperature  Profiles  for  n=.25. 
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Numerical  Solution 


Equation  (1)  was  solved  numerically  after  dividing  the 
material  into  ten  internal  sections  of  equal  width  Ax  as 
shown  in  Figure  7  where  the  labeled  nodes  correspond  to 
specific  locations  within  the  solid. 


Next,  the  finite-difference  equations  were  derived  in  similar 
fashion  as  discussed  on  pages  2-4.  Referring  to  Figure  8, 
the  Taylor  series  expansion  about  Ug  in  region  a  may  be 
found  as  was  (1)>  and  then  rearranged  as 

(AxJ2  3  2Un 


Figure  8.  Node  Arrangement  for  Deriving 
Finite  Difference  Equations 


The  time  derivative  in  region  a  may  be  written  as 


Evaluating  (10)  in  medium  a  about  node  B 


Substituting  (15)  and  (16)  into  (14)  and  rearranging, 


and  similarly  for  region  b 
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It  remains  to  choose  an  approximation  for  Kg  and  , 

the  conductivities  of  region  a  and  region  b,  respectively. 
Four  recommended  approximations  are  the  arithmetic  mean, 
the  harmonic  mean,  the  integral  value,  and  the  value  at  the 
midpoint  of  the  interval.  These  four  approximations  to  the 
value  of  Ka  and  in  (21)  are  included  in  Table  2. 


Table  2 

Approximation  Formulae  for  Conductivity 


Formula 

Name 

K 

a 

(Region  a) 

Kb 

(Region  b) 

Ref 

Arithmetic 

Mean 

ka  +  kb 
- 2 - 

kb  +  Kc 

2 

(15:297) 

Harmonic 

Mean 

2kakb 

2kbkc 

VKc 

(17:45) 

Mid-Point 

k(a;b) 

K  (~7~) 

(12:25) 

Integral 

Value 


l 


(4Xa>U  KflT 


-1 


(Axb) 


Cr  dX 

l  KlxJ 


vl 


(12:24) 


All  expressions  for  conductivity  equate  when  =  Kg  =  K^,  , 
and  (20)  reduces  to  (7)  as  is  required  when  the  diffusivity 
K/pc  ,  is  equal  to  unity  and  the  node  spacing  and  conductivity 
are  constant. 
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Boundary  Condition  at  Infinity 

Fox  (Ref  8:91-92)  suggests  that  an  effective  treatment 
for  the  boundary  condition  at  infinity,  (11)  U(®,t)  =  0  , 

is  to  ensure  the  last  two  nodes  considered  have  zero  values 
to  the  number  of  significant  figures  retained.  This  should 
accurately  satisfy  the  boundary  condition,  and  this  method 
was  used  in  this  investigation. 

Discontinuity 

The  initial  and  boundary  conditions 

U  (0,t)  =  1 

U  (x,0)  *  0 

give  rise  to  a  discontinuity  at  the  origin  because  the  limit¬ 
ing  value  of  the  initial  temperature  is  unity  as  x  tends  to 
zero,  whereas  the  limiting  value  of  the  boundary  temperature 
is  zero  as  t  tends  to  zero.  This  discontinuity  can  give 
rise  to  a  poor  finite-difference  solution  near  x  equal  zero 
for  small  values  of  time  (Ref  20:58).  A  similar  discontinuity, 
in  an  investigation  by  Martin  (Ref  10:72),  caused  poor 
results  when  untreated.  Martin  eliminated  the  discontinuity 
by  substituting  the  exact  analytical  solution  after  the  initial 
time  step  (Ref  10:73)  and  reported  a  significant  improvement 
in  solution  accuracy.  Smith  (Ref  20:58-59)  suggests  a  method 
for  handling  this  problem  which  involves  a  transformation  of 
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the  independent  variables  from  (x,t)  to  (X,T)  where 


The  result  of  this  transformation  is  an  expansion  of  the  origin 
at  x  =  t  =  0  onto  the  positive  side  of  the  new  X  axis, 
while  the  old  x  axis  is  concentrated  at  a  point  at  infinity 
on  the  new  X  axis.  This  transforms  the  jump  condition  at 
x  =  t  =  0  into  a  continuous  change  defined  over  the  positive 
side  of  the  X  axis.  Another  approach,  suggested  by  Ames 
(Ref  1:224-225),  is  to  use  the  mean  initial  value  of  the 
temperature  at  the  origin  for  the  first  time  step  to  reduce 
sensitivity  of  the  finite  difference  approximation  to  the 
discontinuity.  A  fourth  method  discussed  by  Ames  (Ref  1:228) 
and  Ketter  (Ref  10:335-339,  371-379),  and  used  by  Campbell, 
Kaplan,  and  Moore  (Ref  2:325-326)  is  that  of  mesh  refinement. 

In  mesh  refinement,  the  effect  of  the  discontinuity  is  dimin¬ 
ished  by  subdividing  the  space/time  mesh  near  the  discontinuity. 


Error  Analysis 


Background.  Solutions  of  partial  differential  equations 
(PDE)  by  methods  of  finite  difference  equations  (FDE)  are 
subject  to  several  types  of  errors  referred  to  in  this  report 
as  truncation  error,  round-off  error,  and  discretization  error. 
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Truncation  error  is  the  difference  between  the  exact  solution 
of  the  FDE  and  the  analytic  solution  of  the  PDE.  The  truncation 
error  is  due  to  the  use  of  finite  increments  in  approximating 
the  derivatives  and  is  generated  by  the  truncation  of  higher 
order  terms  in  the  approximation  of  the  PDE.  The  order  of  the 
truncation  error  is  approximated  by  the  estimate  of  the  magni¬ 
tude  of  the  first  neglected  term  in  the  Taylor  series  expansion 
used  to  derive  the  FDE  (Ref  1:24).  Round-off  error  is  due  to 
the  limitation  on  the  number  of  significant  figures  carried  by 
the  computer  in  the  solution  of  the  problem.  The  magnitude  of 
this  error  is  a  function  of  the  number  and  type  of  calculations 
involved  in  the  problem  solution.  If  an  infinite  number  of 
decimal  places  were  carried  in  the  computations,  the  round-off 
error  would  vanish  (Ref  19:298).  In  this  report,  the  algebraic 
sum  of  the  truncation  error  and  round-off  error  is  called  the 
discretization  error.  This  error  is  the  difference  between 
the  exact  analytic  solution  and  the  computed  finite  difference 
solution. 

General.  The  accuracy  attainable  with  finite  difference 
methods  will  generally  depend  upon  the  fineness  of  the  grid 
spacing  and  the  order  of  the  finite  difference  approximation. 
Reducing  the  grid  spacing  increases  the  number  of  equations 
to  be  solved,  increasing  the  required  number  of  calculations 
and  the  problem  of  round-off  error.  Using  higher  order  approx¬ 
imations  yields  greater  accuracy  for  a  fixed  mesh  size,  but 
results  in  complicated  finite-difference  expressions  (Ref  6:377). 


20 


Although  it  is  the  discretization  error  that  is 
actually  being  evaluated,  it  is  reasonable  to  assume  this 
error  roughly  equal  to  the  truncation  error.  Therefore, 
if  the  error  is  given  by 

e  =  g ( Ax ) 2 

for  a  node  spacing  of'  Ax  ,  then  the  result  of  decreasing 
the  spacing  to  Ax/2  can  be  expressed  by 

e.  g  (Ax)2 

=  — -  -  4 

h  gjj  (Ax/2)  2 

Similarly,  the  effect  of  halving  the  node  spacing  when  the 
error  is  given  by 

e  =  g  ( Ax  )  4 

is 

e,  g-i  (Ax)  4 

=  Ji -  =  16 

g^  (Ax/2) 4 

The  value  of  the  ratio  of  the  errors  due  to  a  change  in  the 
node  spacing  by  a  factor  of  two  is  defined  as  the  discreti¬ 
zation  error  ratio  in  this  investigation.  Thus,  a  discreti 
zation  error  ratio  of  4  indicates  accuracy  of  order 
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(Ax)2  (O(Ax)2)  ,  and  a  ratio  of  16  indicates  0(Ax)4 
accuracy.  The  truncation  error  associated  with  any  particu¬ 
lar  finite  difference  scheme  may  usually  be  determined  by  a 
Taylor  series  expansion  of  the  finite-differences  about  a 
particular  point  within  the  mesh.  Using  this  method,  Campbell 
(Ref  13:20-22)  has  shown  that,  for  solutions  to  problems 
involving  homogeneous  materials,  the  Explicit,  Pure-Implicit, 
and  Crank-Nicolson  F.D.M.s  are  second  order  accurate,  0(Ax)2  , 
and  the  Douglas  method  is  fourth  order  accurate,  O(Ax)1'  . 

Mitchell  (Ref  12:28)  derives  the  truncation  error  for 
an  implicit  finite-difference  scheme  corresponding  to  the 
Crank-Nicolson  formulation  for  variable  conductivity  and  shows 
it  to  be  0(Ax)2  .  However,  Smith  (Ref  20:216),  in  treating 
irregular  boundaries,  suggests  that  unequal  node  spacing  for 
problems  involving  constant  conductivity  increases  the  trunca¬ 
tion  error  to  O(Ax)  .  The  author  did  not  derive  the  theoreti¬ 
cal  truncation  error  for  finite-differences  involving  variable 
conductivity  and  unequal  node  spacing,  but  suggests  that  it 
would  be  0 (Ax)  on  the  basis  of  the  cited  report  by  Smith. 

Absolute  Error  Magnitude.  The  basic  objective  of  this 
thesis  was  to  compare  the  finite  difference  solutions  obtained 
using  the  Explicit,  Pure-Implicit,  Crank-Nicolson  and  Douglas 
methods.  For  this  analysis,  absolute  error  magnitudes  were 
used.  That  is, 

e  =  UA.S.  "  UF.D.S 
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a-firSTiS 


where 


U 


A.  S. 


U 


F.D.S. 


=  absolute  error  magnitude  which  is  the 

absolute  value  of  the  discretization  error 
=  exact  value  of  U  at  node  x  at 
time  t  from  the  analytic  solution 
=  finite  difference  approximation  of  U 
at  node  x  at  time  t 


Thus,  the  results  of  this  investigation,  as  discussed  in 
Chapter  IV  and  depicted  graphically  in  Appendix  C,  are 
presented  in  terms  of  the  absolute  value  of  the  discreti¬ 
zation  errors  and  the  discret izat ion  error  ratios  for  each 
finite-difference  method  evaluated. 


Stability  Analysis 

Instability.  If  the  magnitude  of  the  difference 
between  the  exact  numerical  solution  and  the  finite- 
difference  approximation  grows  exponentially  as  the  calcu¬ 
lation  proceeds,  then  the  numerical  scheme  is  termed  unstable 

Condition  for  Stability.  For  homogeneous  materials 
where  the  conductivity  is  constant,  the  stability  limits 
are  shown  in  Table  3. 

When  conductivity  varies  only  with  the  spatial  variable, 
as  in  this  investigation,  Richtmyer  and  Morton  (Ref  18:195) 
give  the  general  condition  for  stability  as 
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Table  3 

Stability  Criteria  for  Finite  Difference 

|  Methods 

for  Homogeneous 

Material 

_ -  _ _ _ _ 

Scheme 

f 

Stability  Criteria 

Explicit 

0 

At/ (Ax) 2  £  h 

Pure  Implicit 

1 

Unconditional 

Crank-Nicolson 

h 

Unconditional 

Douglas 

uj  _  (Ax)2,\ 
6At  } 

Unconditional 

General  Equation 

Js  <  f  s  1 

0  s  f  <  h 

Unconditional 

2At/(Ax)zsl/(l-2f) 

2a (x) At  (1  -  2f)/(Ax)2  <  1 


(22) 


where 

a(x)  =  K/pc,  thermal  diffusivity 

and  all  other  variables  have  previously  been  defined.  Using 
Eq  (22)  stability  criteria  may  be  found  for  each  method  by 
appropriate  substitution  for  the  variable  f  as  follows: 
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Explicit  Method: 


Substituting  f  =  0  into  (22) 

2a (x) At/  (Ax) 2  <  1 
and  the  stability  condition  is 
a(x)At/(Ax)2  <  h 

Implicit  Method: 

Substituting  f  =  1  into  (22) 

2a(x)At(l-2)/(Ax)2  <  1 

Thus,  the  requirement  for  stability  is 

-2a (x) At/ (Ax) 2  <  1 

which  is  always  true  as  a(x)  ,  At  and  (ax)  are 
always  positive,  thus  the  Implicit  Method  is 
unconditionally  stable. 

Crank-Nicolson  Method: 

Substituting  f  =  h  into  (22) 

2a (x) At (1-1)/ (Ax) 2  <  1 

Thus  the  requirement  for  stability  is 

0  <  1 

which  is  always  true,  and  the  Crank-Nicolson  Method 
is  unconditionally  stable. 


(23) 
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Douglas  Method: 

Substituting  f  =  %(1  -  — - — )  into  (22) 

6a(x)At 

2a (x) At  1-2(35(1  =  /C^x)2  <  1 

Rearranging 

2a(x)At[  (Ax)2/6a(x)At]/(Ax)2  <  1 
and  the  requirement  for  stability  is 


which  is  always  true.  Thus,  the  Douglas  Method 
is  unconditionally  stable. 

As  was  true  for  homogeneous  material  problems,  only 
the  Explicit  FDE,  for  the  methods  investigated,  presents  a 
problem  concerning  stability.  Substitute 

a (x)  =  K(x)/pc  =  KQxn/pc 

into  (23)  and  let  KQ/pc  =  1.0  .  Then,  Eq  (23)  becomes 

xnAt/(Ax)2  <  H 

which  is  the  stability  limit  for  the  Explicit  FDE  as  used 
in  this  investigation.  For  each  value  of  n  ,  At  ,  and 
Ax  ,  a  maximum  value  of  x  is  calculated  for  which  the 
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Figure  9.  Stability  Curves  for  Explicit 
Finite  Difference  Method  for 
n=0,  . 1,  .25,  .5,  and  .75 


Table  4 


Maximum  Value  of  x  for  Stability 
of  Explicit  Finite  Difference  Method 
where  Criteria  for  Stability  is  xnAt/(Ax)2  <  h 


t/C  x) 


.050 

.100 

.150 

.200 

.250 

.300 

.350 

.400 

.450 

.500 

.550 

.600 

.650 

.700 

.750 

.800 

.850 

.900 

.950 

1.000 


n  =  .  5 


10000. 000 

100.000 

21.544 

625.000 

25.000 

8.550 

169350.878 

123.457 

11.111 

4.979 

9536.743 

39.063 

6.250 

3.393 

1024.000 

16.000 

4.000 

2.520 

165.382 

7.716 

2.778 

1.976 

35.401 

4.165 

2.041 

1.609 

9.313 

2.441 

1.563 

1.347 

2.868 

1.524 

1.235 

1.151 

1.000 

1.000 

1.000 

1.000 

.386 

.683 

.826 

.  881 

.162 

.482 

.694 

.784 

.073 

.350 

.592 

.705 

.035 

.  260 

.  510 

.639 

.  017 

.198 

.444 

.582 

.009 

.153 

.391 

.  534 

.005 

.  120 

.346 

.493 

.003 

.095 

.309 

.457 

.002 

.  077 

.  277 

.425 

.  001 

.  062 

.  250 

.397 
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III.  Procedure 


General  Approach 

Initial  Phase.  Phase  One  of  this  investigation  began 
with  a  study  of  conduction  heat  transfer  and  a  literature 
search  which  included  a  search  for  one  and  two  dimensional 
transient  heat  conduction  problems  involving  variable  con¬ 
ductivity  with  analytic  solutions.  Primary  texts  for  the 
heat  transfer  study  were  Myers  (Ref  IS)  and  Patankar  (Ref  16). 
No  two-dimensional  problems  were  found,  but  one  of  several 
one-dimensional  problems  found  in  Carslaw  and  Jaeger  (Ref  4) 
was  chosen  for  investigation  and  is  presented  in  Chapter  II. 

Second  Phase.  This  phase  consisted  of  deriving  the 
finite  difference  equations  for  non-homogeneous  material 
and  variable  node  spacing  using  engineering  notation,  apply¬ 
ing  the  initial  and  boundary  conditions  of  the  problem,  and 
developing  the  matrix  equations  for  the  numerical  solution 
of  the  problem.  The  matrix  equations  are  presented  in 
Appendix  A. 

Third  Phase.  Phase  Three  began  with  development  of  a 
computer  program  to  numerically  evaluate  the  analytic  sol¬ 
ution  to  the  problem  discussed  in  Chapter  II.  Concurrently, 
the  finite-difference  equations  were  programmed  for  constant 
node  spacing,  and  methods  for  analyzing  and  presenting  data 
were  evaluated.  The  method  of  error  analysis  used  by 
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Campbell,  Kaplan,  and  Moore  (Ref  2)  was  chosen  and  is 
explained  in  Chapter  IV.  Thus,  the  results  are  presented 
as  the  absolute  value  of  the  pointwise  discretization  error 
(at  x  =  .05  and  x  =  0.1)  and  the  discretization  error 
ratio  at  the  same  points  for  each  numerical  method  for 
several  choices  of  time  and  space  increments  as  was  discuss 
in  Chapter  II. 

Fourth  Phase.  The  fourth  phase  of  the  investigation 
began  with  a  study  of  error  and  stability  in  the  finite 
difference  methods.  Richtmyer  and  Morton  (Ref  18)  and  Fox 
(Ref  8)  were  the  primary  texts  used  in  the  study.  The 
stability  limits  as  discussed  in  Chapter  II  were  used  to 
determine  the  time  and  space  parameters  for  the  computer 
runs.  Also  during  this  phase,  a  method  for  approximating 
the  conductivity  as  discussed  in  Chapter  II  was  chosen,  and 

the  computer  codes  written  in  Phase  Three  were  validated. 

\ 

Fifth  Phase.  The  initial  mesh  of  10  nodes  was  sub¬ 
divided  into  20,  40,  and  then  80  nodes.  Correspondingly, 

At  was  reduced  by  factors  of  four  to  keep  the  Fourier 
Modulus  constant.  This  was  done  for  n  equal  to  .25  and 
.5  .  Table  5  shows  the  relationship  of  absolute  time  to 
incremental'  time  for  spatial  subdivisions  considered.  Next 
the  first  cell  was  subdivided  into  two  equal  segments  to 
reduce  the  effect  of  the  discontinuity  at  the  origin  at 
x  =  t  =  0  ,  as  discussed  in  Chapter  II.  This  was  done  for 


31 


Table  S 


to  Incremental  Time 

for  Spatial 

Subdivisions 

Considered 

L  .  _  _ 

Subdivisions 

Real  Time 

20 

Time  Increment 

40 

80 

.01 

1 

4 

16 

.02 

2 

8 

32 

.05 

5 

20 

80 

.1 

10 

40 

160 

.2 

20 

80 

320 

.5 

50 

200 

800 

20,  40,  and  80  nodes,  for  n  =  .25  and  .5  .  The  resulting 
changes  in  degree  of  accuracy  of  each  method  was  analyzed 
and  compared  against  the  previous  results  from  equal  grid 
spacing  for  all  four  finite-difference  methods  under  inves¬ 
tigation. 

The  most  significant  effect  of  the  boundary-condition- 
at-infinity  (11)  is  to  fix  the  number  of  time  steps  a  calcu¬ 
lation  can  proceed,  based  on  the  size  of  the  time  increment, 

At  ,  to  avoid  violating  this  boundary  condition.  In  order 
to  take  larger  time  steps  and  keep  the  same  number  of  nodes, 
it  is  necessary  to  increase  the  size  of  the  spatial  increment, 
Ax  ,  by  an  amount  adequate  to  ensure  the  temperature  at  the 
last  two  nodes  remains  equal  to  zero.  Thus,  the  initial  grid 
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size  with  20  nodes  was  doubled,  the  first  space  increment 
subdivided  into  two  equal  segments,  and  then  subsequent 
runs  were  made  for  n  =  .25  and  .5  for  20,  40,  and  80 
nodes  for  the  Pure-Implicit,  Crank-Nicolson  and  Douglas 
methods  for  At/ (Ax) 2  equal  to  1.0,  2.0  and  5.0.  These 
results  were  also  analyzed  and  compared,  as  were  the  previous 
results. 

Computer  System  and  Programs 

Computer  System.  The  computer  system  used  for  this 
project  is  one  designed  by  the  Control  Data  Corporation. 

It  consists  of  input  and  output  devices,  peripheral  processors, 
and  two  central  processors.  The  central  processors  are  a 
CDC  6613  and  a  CDC  Cyber  74  which  operate  in  parallel.  Each 
has  131,000  60-bit  words  of  central  memory.  Magnetic  disc 
and  drum  storage  were  used  as  temporary  storage  devices. 

Computer  Programs.  Four  major  programs  were  written 
during  this  investigation.  There  were  also  numerous  other 
programs  written  for  error  analysis  and  plotting  purposes. 

The  language  used  for  all  programs  was  FORTRAN  5. 

Although  cost  and  computer  run  time  are  important  in 
most  computer  work,  the  thrust  of  this  investigation  was  in 
the  direction  of  comparing  the  finite-difference  methods 
investigated  for  accuracy  and  stability.  The  programs  which 
were  written  were  designed  for  generality  and  to  make  use 
of  mass  data  storage  and  temporary  files.  Thus,  no  effort 

33 


was  made  to  compare  cost  and  run  times  for  the  various 
options  used.  However,  with  the  exception  of  the  Explicit 
finite-difference  scheme,  all  the  methods  considered  were 
implicit  involving  tri-diagonal  matrices  and  thus  should 
have  been  roughly  comparable  with  respect  to  computation 
time.  Also,  the  Douglas  Method  was  the  most  complicated 
method  to  code  as  the  weighting  factor,  f  ,  had  to  be 
calculated  for  each  node. 


IV.  Results 


Exact  Solution 

As  discussed  in  Chapter  II,  the  analytic  solution, 

Eq  (12),  contained  an  integral  which  had  to  be  evaluated 
numerically.  To  validate  the  computer  program  for  evaluating 
the  analytic  solution,  the  parameter  n  was  set  equal  to 
zero.  Since  a  zero  value  of  n  is  equivalent  to  constant 
conductivity  (since  K  =  KqX11  in  this  investigation),  the 
author  sought  and  found  another  form  of  the  solution  for 
the  investigated  problem  for  constant  conductivity: 


U(x,t)  =  1-erf  ( — - — ) 

/No¬ 


where 

erf  =  error  function 

K  =  conductivity  =  constant 

and  x  and  t  have  been  previously  defined.  An  ASD 
library  routine  was  then  used  to  evaluate  the  error  function. 
A  comparison  of  the  results  for  several  values  of  x  and 
one  value  of  t  is  provided  in  Table  6  as  an  example  of 
the  results  of  this  program  validation. 

Numerical  Solutions 

Conte  and  DeBoor  (Ref  6)  was  the  source  for  the  Thomas 
Algorithm  as  programmed  for  the  tri-diagonal  matrix  solutions 
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A 

Table 

Comparison  Table 

6 

for  Validating 

the  Computer 

Program  for  the 

Exact 

Analytical  Solution 

(t  =  .00625) 

X 

Thesis  Solution 

Library  Solution 

.  05 

. 65472084 

.65472085 

.1 

.37109335 

.37109337 

.  2 

. 07363830 

. 07363827 

.4 

. 00034672 

.00034662 

for  this  investigation.  After  all  programs  were  debugged 
and  error-free,  each  was  modified  to  solve  an  applicable 
heat  transfer  problem  from  Carnahan,  et  al_  (Ref  3)  and  the 
results  were  verified  with  the  published  textbook  results 
to  complete  program  validation. 

Conductivity  Approximation 

The  results  of  the  error  comparisons  made  using  three 
common  methods  for  approximating  the  conductivity  are 
presented  graphically  in  Appendix  B  for  the  Explicit,  Pure- 
Implicit,  Crank-Nicolson,  and  Douglas  finite  difference 
methods  for  a  Fourier  Modulus  of  h  and  with  n  =  .25 
and  n  =  . 5  . 

For  conductivity  which  varies  linearly,  the  mid-point 
and  arithmetic  mean  approximations  are  identical.  For  very 
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small  intervals  and  slowly  varying  conductivity  (non-linear 
variation),  the  mid-point  and  arithmetic  mean  approximations 
should  be  nearly  equal.  When  the  conductivity  varies  by 
relatively  large  amounts,  the  harmonic  mean  should  provide 
the  best  approximation  of  the  three  methods  examined. 

While  the  arithmetic  and  harmonic  mean  approximations 
were  generally  more  accurate  than  the  mid-point  approximation 
for  both  variations  on  the  problem  presented,  (n  =  .25  and 
n  =  .5),  there  was  no  consistency  as  to  whether  the  harmonic 
or  arithmetic  mean  was  the  more  accurate.  Thus,  the  arith¬ 
metic  mean  was  chosen  for  the  rest  of  the  investigation  as 
it  was  the  simpler  method  to  program. 

Discontinuity 

The  discontinuity  at  x  *  t  =  0  as  discussed  in 
Chapter  II  was  treated  by  mesh  refinement  as  discussed  in 
Chapter  III.  The  order  of  accuracy  improved  by  a  factor  of 
two  over  the  unmodified  solutions  for  all  finite-difference 
methods  with  n  equal  to  .25  as  shown  in  Figures  10,  11, 
12,  and  13.  The  improvement  is  not  dramatic.  In  fact,  for 
the  variation  where  n  ,  the  exponent  in  the  conductivity 
term  was  equal  to  0.5,  the  error  improvement  in  the  graded 
mesh  compared  to  the  equally  spaced  mesh  disappeared  after 
approximately  ten  time  steps  which  is  shown  by  comparing 
Figures  14  and  15,  and  16  and  17.  More  improvement  in  error 
reduction  could  be  made,  perhaps  with  one  of  the  other 
methods  discussed  in  Chapter  II. 
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FIGURE  11.  Absolute  value  of  discretization  error  vs 
Time  at  x=.05:  Graded  mesh  n=25 
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FIGURE  12.  Absolute  value  of  discretization  error  vs 
X  at  Time=05.  Ungraded  mesh  n=.25 
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FIGURE  14.  Absolute  value  of  discretization  error  vs 
Time  at  x=1:  Ungraded  mesh  n=5 
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FIGURE  17  Absolute  value  of  discretization  error  vs 
X  at  Ti me=  06:  Graded  mesh  n=.5 
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Although  the  Pure-Implicit,  Crank-Nicolson,  and 
Douglas  Methods  are  unconditionally  stable,  the  graphic 
results  in  Appendix  C,  Section  II,  show  Crank-Nicolson  and 
Douglas  Methods  are  not  unconditionally  accurate.  Patankar 
and  Beliga  (Ref  17:27-31)  and  Crandall  (Ref  7:389)  discuss 
conditions  and  limitations  for  oscillatory  behavior  of 
finite  difference  methods  in  detail.  While  Crandall's 
results  do  not  directly  apply  to  problems  involving  variable 
conductivity,  the  results  of  Appendix  C  imply  that  analogous 
oscillatory  limits  must  exist.  Patankar  and  Beliga  offer 
two  alternative  finite-difference  schemes  which  do  apply 
to  problems  involving  variable  conductivity  and  each  is 
unconditionally  stable  and  has  no  oscillatory  behavior. 
However,  time  and  computer  resources  did  not  permit  using 
these  methods  in  this  investigation.  The  author's  results 
did  verify  that  the  Pure-Implicit  scheme  is  free  from 
oscillatory  behavior  as  was  shown  by  Patankar  and  Beliga 
(Ref  17:29)  . 

Error 

With  the  exception  of  the  results  for  n  =  .5  for  an 
ungraded  mesh,  the  discretization  error  ratio  ranged  from 
1.3  to  1.7  for  all  methods  on  all  variations  of  the  problem. 
This  indicates  0(Ax)  accuracy  as  discussed  in  Chapter  II. 
In  all  graphs  in  Appendix  C,  as  well  as  Figures  10-17,  it 


>r'i*aiai>.gifaitti 


is  noted  that  the  absolute  value  of  the  discretization  error 
and  also  the  discretization  error  ratios  exhibit  erratic 
behavior  at  very  early  times.  This  behavior  with  respect 
to  time  is  characteristic  of  finite-difference  solutions  of 
parabolic  equations,  and  all  plots  show,  as  Smith  (Ref  20:58) 
indicates,  a  characteristic  improvement  in  accuracy  near  x 
equal  zero  as  the  time  increases. 

While  the  graphic  results  are  discussed  in  more  detail 
in  Appendix  C,  a  general  result  was  that  for  n  equal  to 
.25,  the  Douglas  Method  was  usually  the  most  accurate  method 
and  for  n  equal  to  .5,  the  Pure-Implicit  Method  was  the 
most  accurate.  However,  there  was  usually  small  difference 
in  the  relative  accuracy  of  the  methods  with  the  exception 
of  when  the  Crank-Nicolson  and  Douglas  Method  solutions 
exhibited  oscillatory  behavior.  Then,  the  Pure-Implicit 
Method  was  significantly  most  accurate. 

Error  versus  x.  The  discretization  error  ratio  and 
absolute  value  of  the  discretization  error  were  each  plotted 
against  x  for  each  variation  of  the  parameters  n  and 
At/ (Ax) 2  .  Representative  graphs  are  included  in  Appendix  C. 
As  at  x  equal  to  0.05  and  0.1,  the  discretization  error 
ratio  improved  with  time  to  a  range  of  1.3  to  1.7  for  all 
variations  of  the  problem.  All  finite-difference  methods 
were  approximately  equal  in  accuracy  except  when  At/ (Ax) 2 
was  equal  to  5.  Then,  the  Pur e- Implicit  Method  was  signi¬ 
ficantly  most  accurate  across  the  mesh.  Also,  the  solution 
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at  nodes  where  the  temperature  was  very  nearly  zero  behaved 
similarly  to  the  solution  near  x  =  t  =  0  .  That  is,  the 
solution  improved  with  time  at  all  points. 


V.  Conclusions  and  Recommendations 


Conclusions 

Although  the  Crank-Nicolson  and  Douglas  finite- 
difference  schemes  are  increasingly  more  sophisticated  than 
the  traditional  Explicit  and  Pure-Implicit  schemes,  they 
did  not  prove  any  significant  advantage  in  accuracy  in  this 
investigation.  In  all  cases,  for  both  variations  of  the 
problem  (n  =  .25,  n  =  .5)  there  was  little  difference  in 
accuracy  among  the  methods.  Because  of  the  oscillatory 
behavior  of  the  Crank-Nicolson  and  Douglas  schemes  when  the 
value  of  the  Fourier  Modulus  was  large,  the  Pure-Implicit 
scheme  proved  most  advantageous. 

Although  not  verified  theoretically,  all  four  finite 
difference  methods  appear  to  be  first  order  accurate  for 
the  type  of  problem  investigated,  when  the  discontinuity 
at  x  =  t  =  0  is  treated  by  mesh  refinement. 

Recommendat ions 

Martin  (Ref  11)  achieved  a  significant  improvement 
in  the  discretization  error  ratio  when  he  treated  a  dis¬ 
continuity  similar  to  the  one  handled  in  this  investigation 
by  substituting  the  exact  analytic  solution  after  the 
first  time  step.  The  same  should  be  done  for  this  problem 
to  see 'if  similar  improvements  are  achieved.  However,  this 
method  of  treating  the  discontinuity  has  no  real  significant 
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advantage  since  the  analytic  solution  is  rarely  known  for 
"real -world"  engineering  problems.  Similarly,  the  boundary- 
condition-at -inf inity  may  be  treated  by  substitution  of  the 
exact  solution  for  the  last  two  nodes  for  all  time  steps. 
This  would  remove  the  limitation  on  the  number  of  time  steps 
a  calculation  may  proceed  as  discussed  in  Chapter  II. 

The  choice  of  the  conductivity  approximation  was  based 
on  a  comparison  of  absolute  discretization  error  at  x  =  0.1 
for  equal  spaced  nodes  only.  The  majority  of  the  results 
obtained  were  from  calculations  involving  a  graded  network. 
As  shown  in  Appendix  C,  the  solutions  of  all  four  finite- 
difference  methods  investigated  behaved  differently,  depend¬ 
ing  upon  the  value  of  n  ,  the  variable  in  the  exponent  in 
the  conductivity  term,  K  =  KQxn  .  Thus,  comparisons  of  the 
finite-difference  methods  should  be  made  using  the  harmonic 
mean  for  the  conductivity  approximation,  at  least  for 
n  =  .5  .  The  harmonic  mean  may  provide  a  more  satisfactory 
approximation  to  the  exact  solution  than  the  arithmetic 
mean,  since  the  conductivity  varies  more  across  a  given 
interval  for  larger  values  of  n  . 

An  analytic  expression  for  the  truncation  error  should 
be  derived  and  verified  for  finite-differences  for  variable 
conductivity  and  mesh  size.  -These  methods  should  also  be 
tested  on  other  problems  involving  variable  conductivity 
where  the  conductivity  expression  is  simpler  (there  are 
several  suitable  solved  problems  in  Carslaw  and  Jaeger 


(Ref  4:323-324),  and  the  results  compared  with  those  of 
this  investigation. 

Although  Crank-Nicolson  and  Douglas  schemes  have 
the  advantage  of  unconditional  stability,  they  both  have 
oscillatory  limits.  An  investigation  should  be  made  of 
the  use  of  the  exponential  scheme  or  power-law  scheme  as 
discussed  by  Patankar  and  Beliga  (Ref  17:27-32)  for  a 
problem  involving  variable  conductivity,  and  the  results 
included  with  those  of  this  investigation. 
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APPENDIX  A 


Development  of  Matrix  Equations 

The  general  form  of  the  finite-difference  equation, 
as  discussed  in  Chapter  II,  is 


Ax  +  Ax, 

c-V-^  (u;  -  V 


-  fCsr  -  UP  -  sr  <ub  - 

a  a 

K  K 

*  O-OCet  (Ua  -  “«)  -  7PT-  tu»  -  ur)D 
a 


B  Axb  B 


(A- 


which,  for  the  Explicit  method  (f  =  0)  may  be  written  as 


Ug  =  2PUa  +  (1  -  2p  -  2q)  UB  +  2qUc 


where 


and 


At  K 

__  a 

P  ~  pcAxaCAxa+Ax“J 


At  K, 


q  ■ 


pcAxb(Axa+Axb) 
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Using  the  boundary  condition  U(0,t)  =  1  ,  the  following 
matrix  relation  results  (for  four  nodes): 

U/  =  AU  +  C 

where 


For  the  Pure-Implicit  method  (f  =  1),  Eq  (A-l)  may  be 
written 

-2pU'  +  (1+  2p  +  2q)Ug  -  2qU'  =  Ug 
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where  p  and  q  are  as  previously  defined.  The  matrix 
relation  (for  four  nodes)  is 


BU'  +  C  =  U 

where 


For  the  Crank-Nicolson  method  (f  =  %),  Eq  (A-l)  may  be 
written  as 

-pU'  +  (l+p+q)U'  -  qU'  =  pUA  +  (l-p-q)Ug  +  qUc 
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with  p  and  q  as  previously  defined,  and  the  matrix 
relation  written  as 


AIT+  C'  =  BU  +  C 


where 


A 
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For  the  Douglas  Method,  f  =  4(1  -  At^  *  an<^  t^ie 

value  of  f  at  each  node  varies  with  variations  in  the 
conductivity  and  the  grid  spacing.  Thus,  Eq  (A-l)  is 
rewritten 

-2pfAU^  +  C l+2f B (p+q) D  Ug  -  2qfcU£ 

=  2p(l-fA)UA  +  C (f g-1) (p+q) 3  UB  +  2q(l-fc)Uc  (A-2) 

where  p  and  q  have  been  previously  defined  and  fA  *  fg  » 
f^,  denote  the  value  of  the  weighting  factor  f  at  nodes  A  , 

B  ,  and  C  ,  respectively.  Thus,  (A-2)  may  be  rewritten  as 

-2rUA  +  (l+2s)Ug  -  2tU£ 

=  2(p-r)UA  +  (s-p-q)Ug  +  2(q-r)Uc 

where 

r  ‘  PfA 

l  -  qfc 
s  =  (p+q)fg 

and  all  other  variables  have  been  previously  defined. 

The  matrix  relation  (for  four  nodes)  is 

AU'  +  C'  =  BU  +  C 
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APPENDIX  B 


Comparison  of  Conductivity  Approximation  Methods 

This  appendix  contains  the  graphical  results  of  the 
comparison  of  three  conductivity  approximation  methods 
which  were  discussed  in  Chapter  II.  The  results  are  plotted 
as  the  absolute  value  of  the  discretization  error  versus 
time  at  x  =  .1  for  the  Explicit,  Pur e- Implicit ,  Crank- 
Nicolson  and  Douglas  finite  difference  methods,  for  n  =  .25 
and  n  =  .5  .  In  all  cases,  the  numerical  solution  was 
obtained  using  a  mesh  of  40  equally  spaced  nodes  with 
At/ (Ax) 2  ,  the  Fourier  Modulus,  equal  to  0.5. 
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FIGURE  1 

IMPLICIT  F.O.  Method.  Error  comparison  at  x=0.1 
(for  n=.25)for  conductivity  approximation  schemes 


CNJ 


FIGURE  2 

EXPLICIT  F.D  Method:  Error  comparison  at  x=0.1 
(for  n=.25)for  conductivity  approximation  schemes 
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FIGURE  3 

Method:  Error  comparison  at 
conductivity  approximation 


FIGURE  4 

DOUGLAS  F.D.  Method:  Error  comparison  at  x=0.1 
(for  n=25)for  conductivity  approximation  schemes 
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FIGURE  5 

IMPLICIT  F.D.  Method.  Error  comparison  at  x=01 
(for  n=.50)for  conductivity  approximation  schemes 


LEGEND 


□  =  AR  I  TH.  MEAN 
0  =  MID— POINT 
A  =  HARMON  1C  MEAN 


0.000 


0.005 


0.010 


0  015 


0.020 


0.025 


TIME 


FIGURE  6 

EXPLICIT  F.D.  Method:  Error  comparison  at  x=0.1 
(for  n=.50)for  conductivity  approximation  schemes 
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FIGURE  7 

CRANK-N I  COLSON  F.D.  Method:  Error  comparison  at 
x=0.1  (for  n=.50)for  conductivity  approximation 
schemes 


FIGURE  8 

DOUGLAS  F.D.  Method:  Error  comparison  at  x=0.1 
(for  n=.50)for  conductivity  approximation  schemes 
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APPENDIX  C 


Computer  Generated  Plots  of  Results 

This  appendix  contains  a  sample  of  the  graphical 
results  of  this  investigation  presented  as  plots  of  the 
absolute  value  of  the  discretization  error  versus  time  and 
node  location  and  the  discretization  error  ratio  versus  time 
and  node  location  for  the  Explicit,  Pure-Implicit,  Crank- 
Nicolson  and  Douglas  finite-difference  schemes.  For  all 
plots,  the  absolute  value  of  the  discretization  error  is 
labeled  absolute  error.  Discretization  error  and  discreti¬ 
zation  error  ratio  have  been  previously  defined. 

The  plots  are  presented  in  two  sections.  Each  section 
is  prefaced  with  a  short  descriptive  note.  However, 
individual  plots  are  not  necessarily  discussed.  The  sections 
are  organized  according  to  the  value  of  the  Fourier  Modulus 
(At/ (Ax)2)  used  in  calculation.  Section  I  includes  the  graph¬ 
ical  results  of  the  Explicit,  Pure- Implicit ,  Crank-Nicolson 
and  Douglas  Methods  with  the  Fourier  Modulus  equal  to  0.5 
and  n  equal  to  .25.  Section  II  includes  the  graphical 
results  of  the  Pure-Implicit,  Crank-Nicolson  and  Douglas 
Methods  with  the  Fourier  Modulus  equal  to  5  and  n  =  .25. 

The  preface  of  each  of  the  sections  includes  a  key  to  the 
organization  of  the  plots  included  in  the  section. 

Although  calculations  were  performed  for  20,  40,  and 
80  nodes,  the  plots  of  absolute  error  versus  time  and 
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absolute  error  versus  x  for  20,  40,  and  80  nodes  were 
so  similar  that  only  those  plots  for  40  nodes  are  included. 
Exceptions  are  in  Section  II  where  the  solution  for  20  nodes 
behaved  differently  than  that  for  40  and  80  nodes.  There, 
the  graph  for  20  nodes  is  also  included. 

As  discussed  in  Chapter  III,  there  was  small  difference 
in  accuracy  among  the  finite-difference  methods  compared 
except  for  the  oscillatory  behavior  observed  with  the  Fourier 
Modulus  equal  to  2  and  5.  Thus,  the  plots  included  in  this 
appendix  are  presented  as  representative  of  all  of  the  data 
which  was  analyzed  in  this  investigation  as  discussed  in 
Chapter  III. 


Section  I 


This  section  includes  the  graphical  results  of  the 
Explicit,  Pure-Implicit,  Crank-Nicolson  and  Douglas  Methods 
with  the  Fourier  Modulus  equal  to  H  and  n  equal  to  .25. 
The  following  key  shows  which  parameter  options  are  included 
in  this  set  of  graphs. 


Table  C-I 

Key  to  the  Plots  in  Section  1 


Figure 

No . 

Number 
of  Nodes 

X 

Time 

Error* 

C-I-l 

40 

.05 

. 00125- 

.  05 

AE 

C-I-2 

40 

.1 

. 00125- 

.  0625 

AE 

C-I-3 

20/40 

.05 

.00125- 

.05 

DER 

C-I-4 

40/80 

.05 

. 00125- 

.05 

DER 

C-I-5 

20/40 

.1 

. 00125- 

.  0625 

DER 

C-I  -6 

40/80 

.1 

. 00125- 

.  0625 

DER 

C-I  -7 

40 

. 05-.4 

.  025 

AE 

C-I-8 

40 

.  05  - .  4 

.  05 

AE 

C-I-9 

20/40 

.  05-.4 

.  025 

DER 

C-I-10 

40/80 

.05-. 4 

.  025 

DER 

C-I-ll 

20/40 

.05-. 4 

.  05 

DER 

C-I-12 

40/80 

.05-. 4 

.05 

DER 

*AE:  Absolute  Error 

DER:  Discretization  Error  Ratio 
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Although  all  methods  exhibit  approximately  the  same  accuracy, 
the  Douglas  Method  is  generally  most  accurate,  providing 
the  least  discretization  error,  and  the  Explicit,  Crank - 
Nicolson,  and  Pure-Implicit  follow  in  descending  order. 
However,  the  order  is  reversed  when  comparing  discretization 
error  ratio,  as  the  Pure- Implicit  Method  error  is  consistently 
most  improved  by  halving  Ax  ,  and  the  Douglas  Method  error 
least  improved.  Similar  results  were  attained  for  n  equal 
to  .5  except  that  the  Pure-Implicit  and  Crank-Nicolson 
Methods  were  then  most  accurate  and  the  Explicit  and  Douglas 
Methods  provided  the  higher  discretization  error  ratios. 

The  graph"  and  discussion  of  results  in  this  section  are 
equally  rtr  -*■  sentative  of  those  for  a  Fourier  Modulus  equal 
to  1 . 
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FIGURE  C— 1—2-  Absolute  value  of  discretization  errorvs 
Time  at  x=1:  n=  .25  40  nodes 
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FIGURE  C— 1-3.  Absolute  value  of  discretization  error  ratio 
vs  Time  at  x=  05  for  20/40  nodes;  n=  25 


71 


DISC.  ERROR  RATIO 
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FIGURE  C— 1—4.  Absolute  value  of  discretization  error  ratio 
vs  Time  at  x=.05  for  40/80  nodes.  n=  25 
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DISC.  ERROR  RATIO 


FIGURE  C- 1—5.  Absolute  value  of  discretization  error  ratio 
vs  Time  at  x=1  for  20/40  nodes:  n=25 
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FIGURE  C— 1—6.  Absolute  value  of  discretization  error  ratio 
vs  Time  at  x=1  for  40/80  nodes:  n~25 
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ABSOLUTE  ERROR 


FIGURE  C— 1—9.  Absolute  value  of  discretization  errorratio 
vs  X  at  Time=025  for  20/40  nodes:  n=.25 
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FIGURE  C—  I— 11.  Absolute  value  of  discretization  error  ratio 
vs  X  at  Ti me=.05  for  20/40  nodes:  n=.25 


Section  II 


This  section  includes  the  graphical  results  of  the 
Pure-Implicit,  Crank-Nicolson,  and  Douglas  Methods  with 
the  Fourier  Modulus  equal  to  5.0  and  n  equal  to  .25. 

The  following  key  shows  which  parameter  options  are  included 
in  this  set  of  graphs. 


Table  C-II 


Key  to  the 

Plots  in  Section  II 

Figure 

Number 

No. 

of  Nodes 

X 

Time 

Error* 

C-II-1 

20 

.05 

.05-1.0 

AE 

C-II -2 

40 

.05 

.05-1.0 

AE 

C-II-3 

20 

.1 

.05-1.0 

AE 

C-II -4 

40 

.1 

. 05-1.0 

AE 

C-II-5 

20/40 

.05 

.05-1.0 

DER 

C-II-6 

40/80 

.  05 

. 05-1.0 

DER 

C-II-7 

20/40 

.1 

.05-1.0 

DER 

C-II-8 

40/80 

.  1 

. 05-1.0 

DER 

C-II-9 

20 

. 05-1.0 

.25 

AE 

C-II-10 

40 

.05-1.0 

.25 

AE 

C-II-11 

20/40 

.05-1.0 

.  25 

DER 

C-II-12 

40/80 

. 05-1.0 

.25 

DER 

C-II-13 

20 

. 05-1.0 

.5 

AE 

C-II-14 

40 

. 05-1.0 

.5 

AE 

C-II -15 

20/40 

. 05-1.0 

.5 

DER 

C-II -16 

40/80 

. 05-1.0 

.5 

DER 

I 

Absolute  Error 

Discretization 

Error  Ratio 

Both  the  Crank-Nicolson  and  Douglas  Method  solutions 
oscillate  with  time.  The  magnitude  of  the  oscillation  is 
greatest  in  the  Douglas  Method  solution.  Oscillations  are 
observed  only  in  the  solutions  obtained  from  calculations 
made  with  a  mesh  of  20  nodes.  For  meshes  of  40  and  80  nodes, 
all  methods  were  of  the  same  accuracy;  the  Douglas  Method 
most  accurate  across  the  mesh,  the  Pure-Implicit  Method  the 
least  accurate.  Similar  results  were  observed  with  the  value 
of  the  Fourier  Modulus  equal  to  2;  however,  the  observed 
oscillations  damped  sooner  for  the  smaller  value  of  the 
Fourier  Modulus.  Results  of  calculations  performed  with 
n  equal  to  .5  paralleled  those  made  with  n  equal  to  .25 
except  that  for  mesh  sizes  of  40  and  80  nodes,  the  Pure- 
Implicit  Method  was  the  most  accurate  and  the  Douglas 
Method  the  least  accurate,  similar  to  results  discussed  in 
Section  I. 
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ERRO 


FIGURE  C— 11—3.  Absolute  value  of  discretization  error  vs 
Time  at  x=  1:  n=  25  20  nodes 
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DISC.  ERROR  RATIO 
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FIGURE  C— 11—6.  Absolute  value  of  discretization  error  ratio 
vs  Time  at  x=05  for  40/80  nodes.  n=25 
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SC.  ERROR  RATIO 


FIGURE  C- 11-9  Absolute  value  of  discretization  error  vs 
X  at  Time=  25;  n=  25  20  nodes 


FIGURE  C— 11—10.  Absolute  value  of  discretization  error  vs 
X  at  Time=  25  n=.25  40  nodes 


FIGURE  C— 11—11.  Absolute  value  of  discretization  error  ratio 
vs  X  at  Time=  25  for  20/40  nodes;  n~25 


FIGURE  C— 11—12.  Absolute  value  of  discretization  error  ratio 
vs  X  at  Time=25  for  40/80  nodes:  n=  25 
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FIGURE  C— 11—13.  Absolute  value  of  discretization  error  vs 
X  at  Time=.5:  n=25  20  nodes 
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FIGURE  C—  il— 14.  Absolute  value  of  discretization  error  vs 
X  at  Time=  5:  n=  25  40  nodes 
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